﻿// NOTE: All references to clockwise motion are viewed from as
// (x, y) => (0, 1) as indicating going up, (1, 0) as indicating going right

module RectangulateOld
    let IssueError msg = 
        printfn "%s" msg

        (*
    let bin = System.Random()
    //let testMatrixTiles = Array2D.init 4 4 (fun i j -> bin.Next(2))
    let mutable testMatrixTiles = Array2D.init 4 4 (fun i j -> 0)
    testMatrixTiles.[0,1] <- 1
    testMatrixTiles.[0,2] <- 1
    testMatrixTiles.[1,0] <- 1
    testMatrixTiles.[1,1] <- 1
    testMatrixTiles.[1,2] <- 1
    testMatrixTiles.[1,3] <- 1
    testMatrixTiles.[2,0] <- 1
    testMatrixTiles.[2,1] <- 1
    testMatrixTiles.[2,2] <- 1
    testMatrixTiles.[2,3] <- 1
    testMatrixTiles.[3,1] <- 1
    testMatrixTiles.[3,2] <- 1
    let testMatrix = Array2D.init 5 5 (fun i j -> (0, 0))
    
    let convexHull = [|(0,1);(0,2);(0,3);
    (1,3);(1,4);
    (2,4);(3,4);(3,3);
    (4,3);(4,2);(4,1);(3,1);(3,0);(2,0);(1,0);(1,1)|]
    let shape = MatrixToVertexList (TileMatrixToVertexMatrix testMatrixTiles) *)
    let testPerimeter = [| (0,2); (0,3); (1, 3); (1,4); (2,4); (3,4); (4,4); (4,3); (4,2); (4,1); (3,1); (3,0); (2,0); (1,0); (0,0); (0,1) |]
    


    //===========================================================================
    // There are two types of data distictions
    // - Whether it is of tile form or vertex form
    // - Whether it is stored as a matrix or as a list
    //===========================================================================
    // Convertex from tile matrix form to vertex list form
    let MatrixToVertexList (shape : int[,]) =
        [| 0 .. shape.GetLength 0 - 1|] |> Array.map (fun i ->
            [| 0 .. shape.GetLength 1 - 1|] |> Array.filter (fun j -> shape.[i, j] = 1)
            |> Array.map (fun j -> (i, j)))
        //|> Array.filter (fun i -> i.GetLength > 0)
        |> Array.collect (fun x -> x)
            
    // Converts from list form to matrix form
    let ListToMatrix (list : (int*int)[]) = 
        let width = list |> Array.map (fun tuple -> fst tuple) |> Array.max
        let height = list |> Array.map (fun tuple -> snd tuple) |> Array.max
        let mutable m = Array2D.zeroCreate (width + 1) (height + 1)
        for i = 0 to list.Length - 1 do
            m.[fst list.[i], snd list.[i]] <- 1
        m
        
    // Converts from tile matrix from to vertex matrix form
    let TileMatrixToVertexMatrix (shape : int[,]) =
        let width = shape.GetLength(0) - 1
        let height = shape.GetLength(1) - 1
        Array2D.init (width + 2) (height + 2) (fun i j ->
            shape.[max (i - 1) 0, max (j - 1) 0] |||
                shape.[min i width, max (j - 1) 0] |||
                shape.[max (i - 1) 0, min j height] |||
                shape.[min i width, min j height]
            ) // 1 if all four corners are 1, 0 otherwise

    // Converts from vertex matrix form back to tile matrix form
    let VertexMatrixToTileMatrix (shape : int[,]) =
        let width = shape.GetLength(0)
        let height = shape.GetLength(1)
        Array2D.init (width - 1) (height - 1) (fun i j ->
            shape.[i, j] * shape.[i + 1, j] * shape.[i, j + 1] * shape.[i + 1, j + 1]
            ) // 1 if all four corners are 1, 0 otherwise
            
    // Prints out a Array2D with x and y coordinates in the expected orientation
    let print2D (ls : int[,]) =
        for y = 0 to (Array2D.length2 ls) - 1 do
            for x = 0 to (Array2D.length1 ls) - 1 do
                printf "%d" ls.[x, y]
            printfn ""

    // Debugging function for printing out a vertex list form
    let printVertexArrayList (list : (int*int)[][]) =
        list |> Array.map (fun i -> printfn "%A\n" <| ListToMatrix i)

    // An implementation of merge sort seems to be unnecessary
    // Kept incase the need for a stable sort is required
    // Needs to be checked to ensure that this implementation is in fact stable
    // Use would be within GrahamScan for the ordering of elements by array, which even now is slightly buggy
    (*
    /// Merge Sort - Sorts ls weighing each element by the function op and performing a less than or equal comparison
    let rec MergeSort op ls =
        match ls with
        | [||] -> [||]
        | [|x|] -> [|x|]
        | ls ->
           // printfn "mergeSort3: %A" ls
            
            // Places the list into a tree data structure
            let rec Split ls n =
                match ls with
                | [||] -> [||],[||]
                | ls when n=0 -> ([||],ls)
                | _ ->
                    let l1, l2 = Split ls.[1 .. ls.Length - 1] (n-1)
                  //  printf "SplitList - %A : %A" l1 l2
                    ((Array.append [|ls.[0]|] l1), l2)
 
            let rec Merge ls1 ls2 op =
                match ls1, ls2 with
                | [||],[||] -> [||]
                | [||],ls2 -> ls2
                | ls1,[||] -> ls1
                | s1, s2 when (op s1.[0]) <= (op s2.[0]) -> Array.append [|s1.[0]|] (Merge s1.[1 .. s1.Length - 1] ls2 op)
                | s1, s2 -> Array.append [|s2.[0]|] (Merge ls1 s2.[1 .. s2.Length - 1] op)

            let ls1, ls2 = Split ls (ls.Length / 2)
            Merge (MergeSort op ls1) (MergeSort op ls2) op *)


    // Performs a graham scan to find the convex hull, includes the colinear points along the perimeter
    // Contains the first point duplicated in the end
    // Contains a single hack such that it is only intended to be used by TracePerimeter
    let GrahamScan (vertexList : (int * int)[]) = 
        // A helper function to swap the two elements index i1 and index i2 in the Array ls
        let swap i1 i2 ls =
            [| 0 .. (Array.length ls) - 1 |]
                |> Array.map (fun i -> match i with | i when i=i1 -> i2 | i when i=i2 -> i1 | _ -> i)
                |> Array.map (fun i -> ls.[i])
        // A helper function that tells whether three points are counter-clockwise or not through the use of the deterimant
        // The three points are a counter-clockwise turn if < 0, clockwise if > 0, and collinear if = 0
        let findOrientation (p1 : int * int) (p2 : int * int) (p3 : int * int) =
           (fst p2 - fst p1) * (snd p3 - snd p1) - (snd p2 - snd p1) * (fst p3 - fst p1)


        // Sorts the vertexList counter clockwise
        let startY =  vertexList |> Array.filter (fun y -> snd y = (vertexList |> Array.map (fun tuple -> snd tuple) |> Array.min))
        let head = startY |> Array.filter (fun x -> fst x = (startY |> Array.map (fun tuple -> fst tuple) |> Array.min))
        let rest =
            //vertexList |> Array.filter (fun x -> x <> head.[0])
            //|> MergeSort (fun i -> (float <| fst i) / (float <| fst i * fst i + snd i * snd i))
            
            [| 0 .. vertexList.Length - 1 |]
            // Removes head from the list
            |> Array.filter (fun i -> vertexList.[i] <> head.[0])
            // Sorts the vertexList by finding the cosine of the angle formed from head to a point with the x-axis
            |> Array.sortBy (fun i ->
                let x = fst vertexList.[i] //- fst head.[0] // Possible fix to the problem that points directly y-positive above head
                let y = snd vertexList.[i] - snd head.[0]   // all cause cosine to squash them together
                (float <| x) / (sqrt (float <| x * x + y * y)))
            |> Array.map (fun i -> vertexList.[i]) 

        let mutable hull = Array.append (Array.append head rest) head
        //printfn "%A" hull
    
        // Begin sorting points starting from the head, making order decision based off the whether clockwise turns are made
        let mutable i = 3
        let mutable M = 2 // M denotes the number of points on the convex hull
        while i < hull.Length do
            // Find next valid point on convex hull
            while (findOrientation hull.[M - 1] hull.[M] hull.[i]) >= 0 do
                // Check if first points are collinear, if so ignore unnecessary points
                if M = 2 then
                    hull <- swap M i hull
                    i <- i + 1
                    //if i >= points.Length then // If it's just a line
                    //    [| |]
                else
                    M <- M - 1

            M <- M + 1
            hull <- swap M i hull
            i <- i + 1
        hull.[0 .. M]

    // Returns in vertex-form list
    // OPTIMIZATION have TracePerimeter also return the characterization information
    // so as to reduce the number of steps done in Rectangulate
    let TracePerimeter (polygon : int[,]) = 
        // Categorizes the vertex vertexFormPoint based on its location in the polyon tileFormMatrix.  If it returns
        //  - 1: A convex point
        //  - 2: A point on the edge
        //  - 3: A concave point
        //  - 4: Point that lies on the interior of the function
        let CategorizePoint (vertexFormPoint : (int * int)) (tileFormMatrix : int[,]) = 
            let x = fst vertexFormPoint
            let y = snd vertexFormPoint
            tileFormMatrix.[x - 1, y - 1] + tileFormMatrix.[x - 1, y] +
                tileFormMatrix.[x, y - 1] + tileFormMatrix.[x, y]
        let add a b = fst a + fst b, snd a + snd b // Adds two vectors in tuple form together
            
        // Creates from the polygon the same polygon matrix except with a padding of 1 square
        // so that cateorize point can categorize effectively
        let paddedWidth = (polygon.GetLength 0) + 1
        let paddedLength = (polygon.GetLength 1) + 1
        let paddedPolygon = Array2D.init (paddedWidth + 1) (paddedLength + 1)  (fun i j ->
            if i = 0 || j = 0 || i = paddedWidth || j = paddedLength then 0
            else polygon.[i - 1, j - 1])
        let vertexList = polygon |> TileMatrixToVertexMatrix |> MatrixToVertexList |> Array.map (fun x -> fst x + 1, snd x + 1)
        //let vertexList = paddedPolygon |> TileMatrixToVertexMatrix |> MatrixToVertexList // forces sequentiality
        


        // From the convex hull of the padded polygon, fill in the perimeter
        // The following commented section is a more sequential way of doing things thorough a for loop over the convexHull array
        let convexHull = GrahamScan vertexList
        // A sequential method for filling it in
        (*let mutable dir = (-1, 0) // for any polygon, beginning from the lowest y, lowest x coordinate, the only possible direction is up
                                // moving in a clockwise fashion
        let mutable perimeter = [| convexHull.[0] |]
        vertexList |> ListToMatrix |> print2D
        
        for i = 1 to convexHull.Length - 1 do
            let mutable cur = perimeter.[perimeter.Length - 1]
            printfn ""
            while cur <> convexHull.[i] do
                //printf "(%d,%d)" (fst cur) (snd cur)
                let t = CategorizePoint cur paddedPolygon
                if t = 1 then   dir <- (snd dir, -fst dir) // Convex turn clockwise
                if t = 3 then   dir <- (-snd dir, fst dir) // Concave turn counter-clockwise
                if t = 4 then   printfn "Error: TracePerimeter has entered into the shape"
                //printf "=>%d %d,%d=>" t (fst dir) (snd dir)
                cur <- add cur dir
                perimeter <- Array.append perimeter [| cur |]
        // Returns the first point repeated as a duplicate in the end
        perimeter.[0 .. perimeter.Length - 1]*)
        
        //printfn "%A" convexHull
        //vertexList |> ListToMatrix |> print2D
      
        // The following code is the same as above however is a more efficient method
        // that break the calculation up into smaller pieces over perhaps several cores
        // Seperate the task in to the smaller tasks of tracing out the space from one point on the hull to the next
        let perimeter =
            convexHull |> Seq.pairwise
            |> Seq.map (fun tuple ->
                // A and B are part of the convex hull and must both be convex
                let a = fst tuple 
                let b = snd tuple
                let x = fst a
                let y = snd a

                // Find the direction to intially undertake
                // Takes the vector from a to b and finds the correct direction based on the tile
                // to which a is the corner point
                // OPTIMIZATION: set this up with a series of logic statements instead
                let generalDir = fst b - x, snd b - y
                let initialDir =
                    if paddedPolygon.[x - 1, y - 1] = 1 && fst generalDir < 0 then (-1, 0)
                    else if paddedPolygon.[x - 1, y - 1] = 1 && snd generalDir < 0 then (0, -1)
                    else if paddedPolygon.[x - 1, y] = 1&& fst generalDir < 0 then (-1, 0)
                    else if paddedPolygon.[x - 1, y] = 1 && snd generalDir > 0 then (0, 1)
                    else if paddedPolygon.[x, y - 1] = 1 && fst generalDir > 0 then (1, 0)
                    else if paddedPolygon.[x, y - 1] = 1 && snd generalDir < 0 then (0, -1)
                    else if paddedPolygon.[x, y] = 1 && fst generalDir > 0 then (1, 0)
                    else (* if tileFormMatrix.[x, y] && snd generalDir > 0 then*) (0, 1)

                let mutable cur = add a initialDir
                let mutable inBetween = [| a |]
                let mutable dir = initialDir
                while cur <> b do
                    let t = CategorizePoint cur paddedPolygon
                    if t = 1 then   dir <- (snd dir, -fst dir) // Convex turn clockwise
                    if t = 3 then   dir <- (-snd dir, fst dir) // Concave turn counter-clockwise
                    if t = 4 then   printfn "Error: TracePerimeter has entered into the shape"
                    //printf "=>%d %d,%d=>" t (fst dir) (snd dir)
                    inBetween <- Array.append inBetween [| cur |]
                    cur <- add cur dir
                inBetween)
            |> Seq.toArray
            |> Array.collect (fun i -> i)
      
        // Now to remove the padding initially introduced padding around the perimeter
        perimeter |> Array.map (fun tuple -> fst tuple - 1, snd tuple - 1)

        
        

    // For testing GrahamScan and TracePerimeter
    (* let test2 = [| (3,0); (15,2); (12,6); (0,4); (12,0); (15,4); (3,6); (0,2); (6,6); (6,4); (9,4); (9,6) |]
    let test = [| (3,0); (9,2); (6,6); (0,4); (6,0); (9,4); (3,6); (0,2) |]
    test |> 
    //let test = [| (3, 0); (0,2); (0,4) |]
    test2 |> ListToMatrix |> print2D
    test2 |> GrahamScan

    
    let mutable testMatrix1 = Array2D.init 4 4 (fun i j -> 0)
    testMatrix1.[0,1] <- 1
    testMatrix1.[0,2] <- 1
    testMatrix1.[1,0] <- 1
    testMatrix1.[1,1] <- 1
    testMatrix1.[1,2] <- 1
    testMatrix1.[1,3] <- 1
    testMatrix1.[2,0] <- 1
    testMatrix1.[2,1] <- 1
    testMatrix1.[2,2] <- 1
    testMatrix1.[2,3] <- 1
    testMatrix1.[3,1] <- 1
    testMatrix1.[3,2] <- 1
    testMatrix1 |> print2D
    
    let mutable testMatrix2 = Array2D.init 4 8 (fun i j -> 0)
    testMatrix2.[0,0] <- 1
    testMatrix2.[0,1] <- 1
    testMatrix2.[0,2] <- 1
    testMatrix2.[0,3] <- 1
    testMatrix2.[0,4] <- 1
    testMatrix2.[0,5] <- 1
    testMatrix2.[0,6] <- 1
    testMatrix2.[0,7] <- 1
    testMatrix2.[1,0] <- 1
    testMatrix2.[1,1] <- 1
    testMatrix2.[1,2] <- 1
    testMatrix2.[1,3] <- 1
    testMatrix2.[1,4] <- 1
    testMatrix2.[1,5] <- 1
    testMatrix2.[1,6] <- 1
    testMatrix2.[1,7] <- 1
    testMatrix2.[2,0] <- 1
    testMatrix2.[2,1] <- 1
    testMatrix2.[3,0] <- 1
    testMatrix2.[3,1] <- 1
    testMatrix2.[2,4] <- 1
    testMatrix2.[2,5] <- 1
    testMatrix2.[3,4] <- 1
    testMatrix2.[3,5] <- 1
    testMatrix2 |> print2D*)


    let SplitPathBetweenPoints a b (list : (int*int)[]) =
        let posA = [| 0 .. list.Length - 1|] |> Array.fold (fun acc elem -> if a = list.[elem] then elem else acc) -1
        let posB = [| 0 .. list.Length - 1|] |> Array.fold (fun acc elem -> if b = list.[elem] then elem else acc) -1
        if posA < posB then
            ((Array.append list.[0 .. posA] list.[posB .. list.Length - 1]), list.[posA .. posB])
        else
            (list.[posB .. posA], (Array.append list.[posA .. list.Length - 1] list.[0 .. posB]))

    // Defines a function that will check if a point is within a list
    let pointIsInList (list : (int*int)[]) (point : (int*int)) =
        (((-1, -1), list) ||> Array.fold (fun acc elem ->
            if point = elem then
                elem
            else
                acc)) <> (-1, -1)

    let pointIsInRectangle (rectangles : ((int*int)[])[]) (point : (int*int)) =
      //  rectangles |> Seq.map (fun i -> pointIsInList point i)
      //      |> Seq.fold (fun acc elem -> acc || elem) false
        let matches = rectangles |> Array.map (fun i -> pointIsInList i point)
        [| 0 .. rectangles.Length - 1|] |> Array.filter (fun i -> matches.[i])

    // Finds the concave vertices from the in-order clockwise outline of the perimeter
    // Does this by tracing points from 
    let FindConcaveVetices (tileShape : int[,]) (perimeter : (int*int)[]) = 
        let vertexShape = MatrixToVertexList <| TileMatrixToVertexMatrix tileShape
        let CorrectBounds a = (a + perimeter.Length) % perimeter.Length
        let AddPoints pt1 pt2 = (fst pt1 + fst pt2, snd pt1 + snd pt2)
        let SubPoints pt1 pt2 = (fst pt1 - fst pt2, snd pt1 - snd pt2)
        [| 0 .. perimeter.Length - 1|] |> Array.filter (fun index ->
            let vertex = perimeter.[index]
            let indexNext = CorrectBounds (index + 1)
            let indexPrevious = CorrectBounds (index - 1)
            if index < 0 then
                IssueError "Vertex not found"

            let fromBack = SubPoints vertex perimeter.[indexPrevious]
            let forward = SubPoints perimeter.[indexNext] vertex

            //printfn "%d<-%d->%d" indexPrevious index indexNext
            //printfn "%A<-%A->%A" convexHull.[indexPrevious] convexHull.[index] convexHull.[indexNext]
            //printfn "%A<->%A" fromBack forward
            (forward <> fromBack) && (pointIsInList vertexShape (AddPoints vertex fromBack)) &&
                (pointIsInList vertexShape (SubPoints vertex forward)))
        |> Array.map (fun i -> perimeter.[i])










    // The following two functions are being worked on to fix a potential bug that I noticed might happened
    // Code is incomplete
    let DoesBetweenExist (perimeter : (int * int)[]) (a : int*int) (b : int*int) =
        let x1 = min (fst a) (fst b)
        let x2 = max (fst a) (fst b)
        let y1 = min (snd a) (snd b)
        let y2 = max (snd a) (snd b)
        (perimeter |> Array.filter (fun p ->
            (p <> a) && (p <> b) && (x1 <= fst p) && (fst p <= x2) && (y1 <= snd p) && (snd p <= y2))).Length > 0

    let IsOnRectangle (perimeter : (int * int)[]) (point : int*int) = 
        let applySearch f1 f2 =
            let list = perimeter |> Array.filter (fun tuple ->
                f1 point = f1 tuple && point <> tuple) |> Array.sortBy (fun tuple -> abs <| f2 point - f2 tuple)
            // Two cases:
            // 1. point is on the perimeter already in which case there is only one element
            // 2. Point is not on the perimeter and there are two points we must consider
            list.[0 .. max 1 (list.Length - 1)] |> Array.filter (fun tuple ->
                not ((point, tuple) ||> DoesBetweenExist perimeter))
        let horizontal = applySearch fst snd
        let vertical = applySearch snd fst
        printfn "%A\n%A" horizontal vertical
        horizontal.Length > 0 && vertical.Length > 0
        









    
    let ChordlessRectangulate (perimeter : (int*int)[]) (shape : int[,]) = 
        let concaveVertices = FindConcaveVetices shape perimeter
        let listOfVertices = [| |]
        //let direction = (fst perimeter.[1] - fst perimeter.[0], snd perimeter.[1] - snd perimeter.[0])
        
        let CorrectBounds a = (a + perimeter.Length) % perimeter.Length
        let AddPoints pt1 pt2 = (fst pt1 + fst pt2, snd pt1 + snd pt2)

        // Gets the direction in which to draw lines
        let taggedCVs = concaveVertices |> Array.map (fun vertex ->
            let index = (-1, [|0 .. perimeter.Length - 1|]) ||> Array.fold (fun acc elem ->
                if vertex = perimeter.[elem] then
                    elem
                else
                    acc)
            //let indexNext = CorrectBounds (index + 1)
            let indexPrevious = CorrectBounds (index - 1)
            if index < 0 then
                IssueError "Vertex not found"

            let fromBack = (fst perimeter.[index] - fst perimeter.[indexPrevious], snd perimeter.[index] - snd perimeter.[indexPrevious])
            //let forward = (fst perimeter.[indexNext] - fst perimeter.[index], snd perimeter.[indexNext] - snd perimeter.[index])
            //let forward = AddPoints (perimeter.[indexNext]) (-1 * fst perimeter.[index], -1 * snd  perimeter.[index])
            
            (fromBack, vertex))

        printf "begin drawing"
        // Draws the lines from convex holes to walls
        let mutable rectangleList = [| perimeter |]
        for i = 0 to taggedCVs.Length - 1 do
            let concaveStart = AddPoints (snd taggedCVs.[i]) (fst taggedCVs.[i])
            let mutable point = concaveStart

            // Draws a cord from the concave vertex to the first point it finds that exists in the rectangle
            // Drawn in the direction from the previous point
            let mutable chord : (int*int)[] = [| |]
            printf "\nentering while loop"
            while (pointIsInRectangle rectangleList point).Length = 0 do
                chord  <- Array.append chord [| point |]
                point <- AddPoints point (fst taggedCVs.[i])
                printf "\n\tchord: %A" chord
                printf "\n\tpoint: %A" point
                printf "\n\trectangleLIst: %A" rectangleList
                printf "\n\tpointIsInRectangle: %A" (pointIsInRectangle rectangleList point)
                printf "\n\tlength: %d" (pointIsInRectangle rectangleList point).Length

            // Searches for the rectangle that contains both the concave vertex we started at,
            // and the point at the end of the chord we have just drawn
            let locA = pointIsInRectangle rectangleList (snd taggedCVs.[i])
            let locB = pointIsInRectangle rectangleList point
            printf "\ndoing rectToSplit"
            let rectToSplit = locA |> Array.filter (fun i ->
                locB |> Array.fold (fun acc elem -> acc || (i = elem)) false)

            // Let's us know if there exists a rectangle to search or not
            if rectToSplit.Length = 0 then
                IssueError "Rectangle with endpoint and concave start not found"
            let index = rectToSplit.[0] // The index of the rectangle to split

            // Remove the rectangle to be split and adds the next to rectangles
            let splitRect = (SplitPathBetweenPoints (snd taggedCVs.[i]) point rectangleList.[index])
            rectangleList <- Array.append rectangleList [| Array.append (fst splitRect) chord |]
            rectangleList <- Array.append rectangleList [| Array.append (snd splitRect) (Array.rev chord) |]
            rectangleList <- Array.append rectangleList.[0 .. index - 1] rectangleList.[index + 1 .. rectangleList.Length - 1]
        rectangleList
        (*
    let testPerimeter2 = [| (0,2); (0,3);
                            (1, 3); (2,3); (3,3); (4,3); (4,4); (4,5);
                            (3,5); (2,5); (1,5);
                            (0,5); (0,6); (0,7);
                            (0,8); (1,8); (2,8); (3,8); (4,8);
                            (5,8); (5,7);
                            (6,7); 
                            (6,6); (7,6); (7,5); (8,5); (8,4); (8,3); (8,2); (7,2);
                            (6,2); 
                            (6,1); (5,1); (4,1);
                            (4,0); (3,0); (2,0); (1,0); (0,0); (0,1) |]
    let concavePerimeter2  = [| (4,3); (4,5); (5,7);(6,6); (7,5); (6,2); (4,1)|]
    printfn "%A" <| ListToMatrix testPerimeter2
    printfn "%A" <| ListToMatrix concavePerimeter2
    *)

    // Rectangulate
    let RectangulateOld (rect : (int*int)[]) (tiledShape : int[,]) =
        let concaveVertices = FindConcaveVetices tiledShape rect
        let applySearch list f0 f1 = 
            let likelyMatchPairs =
                list |> Array.map (fun from ->
                    list |> Array.filter (fun till ->
                        f0 from = f0 till)
                    //|> Array.sortBy (fun tuple -> abs(f1 from - f1 tuple))
                    |> Array.filter (fun till -> // Removes points that 
                        let t = f0 from // Index along the axis shared by both points
                        let a = f1 from // Index along axis on which the points differ
                        let b = f1 till // Index along axis on which the points differ
                        if a < b then // Checks if all the poitns in the two positions
                            [| a + 1 .. b - 1 |] |> Array.fold (fun acc elem ->
                                acc && not (pointIsInList rect (t, elem))) true
                        else
                            [| b + 1 .. a - 1 |] |> Array.fold (fun acc elem ->
                                acc && not (pointIsInList rect (t, elem))) true
                        )
                    |> Array.sortBy (fun tuple -> f1 tuple)) // Organize
                |> Array.filter (fun coordList -> coordList.Length > 1) // Removes the no possilibty matches
            
            [|0 .. likelyMatchPairs.Length - 1|] |> Array.filter (fun i ->
                    [|0 .. i - 1|] |> Array.fold (fun acc elem -> acc && likelyMatchPairs.[i] <> likelyMatchPairs.[elem]) true)
                |> Array.map (fun i -> likelyMatchPairs.[i])
        let horizontal = applySearch concaveVertices fst snd
        let vertical = applySearch concaveVertices snd fst
        
        //printfn "%A\n%A" horizontal vertical
        // Finds the chords between each
        let findMatches (list : ((int * int))[][]) (list2 : ((int * int))[][]) f0 f1 =
            list |> Array.map (fun tupleList ->
                let w = f0 tupleList.[0]
                let a = f1 tupleList.[0]
                let b = f1 tupleList.[1]
                [| 0 .. list2.Length - 1 |] |> Array.filter (fun i ->
                    let u = f1 list2.[i].[0]
                    let c = f0 list2.[i].[0]
                    let d = f0 list2.[i].[1]
                    //if (not ((a < u) && (u < b) && (c < w) && (w < d))) then
                    //printfn "%A<->%A->%b" tupleList list2.[i] ((a < u) && (u < b) && (c < w) && (w < d))

                    ((a < u) && (u < b) && (c < w) && (w < d))
                    ))
        let matches = Array.append (findMatches horizontal vertical fst snd |> Array.map (fun i ->
            i |> Array.map (fun j -> j + horizontal.Length))) (findMatches vertical horizontal snd fst)
        let matchesH = findMatches horizontal vertical fst snd
        //printfn "\n\nbreak\n"
        let matchesV = findMatches vertical horizontal snd fst
            
        //printfn "%A\n\n%A\n\n%A" matches matchesH matchesV
        let independentSet = Rectangulate.MinimumIndependentSet matchesH matchesV
        //let independentSet = MinimumIndependentSet matchesH matchesV
        //printfn "%A" independentSet
        //let independentSet = [|1;2|]
        let chordsToRectangulateOver = independentSet |> Array.map (fun i -> (Array.append horizontal vertical).[i])
        //printfn "%A" (Array.append horizontal vertical)
        
        let mutable rectList : (int * int)[][] = [| rect |]
        for i = 0 to chordsToRectangulateOver.Length - 1 do
            let chord = chordsToRectangulateOver.[i]
            let a = chord.[0]
            let b = chord.[1]
            let dirTemp = (fst b - fst a, snd b - snd a)
            let chordToAdd =
                if fst dirTemp = 0 then
                    [|snd a + 1 .. snd b - 1|] |> Array.map (fun i -> (fst a, i))
                else
                    [|fst a + 1 .. fst b - 1|] |> Array.map (fun i -> (i, snd a))
            //let start = a + dirTemp

            let rList = rectList // Get a immutable snapshot of the current state
            // Find the index of the rectangle that contains both a and b
            let p = (-1, pointIsInRectangle rList a) ||> Array.fold (fun acc elem ->
                if (pointIsInRectangle rList b |> Array.fold (fun acc2 elem2 -> acc2 || (elem2 = elem)) false) then
                    elem
                else
                    acc)
            let splitRect = SplitPathBetweenPoints a b rectList.[p]
            rectList <- Array.append rectList [| (Array.append (fst splitRect) chordToAdd) |]
            rectList <- Array.append rectList [| (Array.append (snd splitRect) (Array.rev chordToAdd)) |]
            //printfn "%A" rectList
            rectList <- rectList.[p + 1 .. rectList.Length - 1]
        rectList |> Array.map (fun chordlessRect ->
            ChordlessRectangulate chordlessRect tiledShape)
        |> Array.collect (fun x -> x)


    // OPTIMIZATION: Have TracePerimeter potentially pass characterization information over so that
    // there's no need to invoke FindConcaveVertices for the Rectangulate routine
    let Rectangulate2 (tiledShape : int[,]) =
        printf "tracing perimeter"
        let rect : (int*int)[] = TracePerimeter tiledShape
        printf "finding concave vertices"
        let concaveVertices = FindConcaveVetices tiledShape rect
        let applySearch list f0 f1 = 
            let likelyMatchPairs =
                list |> Array.map (fun from ->
                    list |> Array.filter (fun till ->
                        f0 from = f0 till)
                    //|> Array.sortBy (fun tuple -> abs(f1 from - f1 tuple))
                    |> Array.filter (fun till -> // Removes points that 
                        let t = f0 from // Index along the axis shared by both points
                        let a = f1 from // Index along axis on which the points differ
                        let b = f1 till // Index along axis on which the points differ
                        if a < b then // Checks if all the poitns in the two positions
                            [| a + 1 .. b - 1 |] |> Array.fold (fun acc elem ->
                                acc && not (pointIsInList rect (t, elem))) true
                        else
                            [| b + 1 .. a - 1 |] |> Array.fold (fun acc elem ->
                                acc && not (pointIsInList rect (t, elem))) true
                        )
                    |> Array.sortBy (fun tuple -> f1 tuple)) // Organize
                |> Array.filter (fun coordList -> coordList.Length > 1) // Removes the no possilibty matches
            
            [|0 .. likelyMatchPairs.Length - 1|] |> Array.filter (fun i ->
                    [|0 .. i - 1|] |> Array.fold (fun acc elem -> acc && likelyMatchPairs.[i] <> likelyMatchPairs.[elem]) true)
                |> Array.map (fun i -> likelyMatchPairs.[i])
        printf "applying search horizontally"
        let horizontal = applySearch concaveVertices fst snd
        printf "apply search vertically"
        let vertical = applySearch concaveVertices snd fst
        
        //printfn "%A\n%A" horizontal vertical
        // Finds the chords between each
        let findMatches (list : ((int * int))[][]) (list2 : ((int * int))[][]) f0 f1 =
            list |> Array.map (fun tupleList ->
                let w = f0 tupleList.[0]
                let a = f1 tupleList.[0]
                let b = f1 tupleList.[1]
                [| 0 .. list2.Length - 1 |] |> Array.filter (fun i ->
                    let u = f1 list2.[i].[0]
                    let c = f0 list2.[i].[0]
                    let d = f0 list2.[i].[1]
                    //if (not ((a < u) && (u < b) && (c < w) && (w < d))) then
                    //printfn "%A<->%A->%b" tupleList list2.[i] ((a < u) && (u < b) && (c < w) && (w < d))

                    ((a < u) && (u < b) && (c < w) && (w < d))
                    ))
        printf "finding matches"
        let matches = Array.append (findMatches horizontal vertical fst snd |> Array.map (fun i ->
            i |> Array.map (fun j -> j + horizontal.Length))) (findMatches vertical horizontal snd fst)
        let matchesH = findMatches horizontal vertical fst snd
        //printfn "\n\nbreak\n"
        let matchesV = findMatches vertical horizontal snd fst
            
        //printfn "%A\n\n%A\n\n%A" matches matchesH matchesV
        printf "finding independentset"
        let independentSet = Rectangulate.MinimumIndependentSet matchesH matchesV
        //let independentSet = MinimumIndependentSet matchesH matchesV
        //printfn "%A" independentSet
        //let independentSet = [|1;2|]
        printf "another IS mapping"
        let chordsToRectangulateOver = independentSet |> Array.map (fun i -> (Array.append horizontal vertical).[i])
        //printfn "%A" (Array.append horizontal vertical)
        
        let mutable rectList : (int * int)[][] = [| rect |]
        for i = 0 to chordsToRectangulateOver.Length - 1 do
            let chord = chordsToRectangulateOver.[i]
            let a = chord.[0]
            let b = chord.[1]
            let dirTemp = (fst b - fst a, snd b - snd a)
            let chordToAdd =
                if fst dirTemp = 0 then
                    [|snd a + 1 .. snd b - 1|] |> Array.map (fun i -> (fst a, i))
                else
                    [|fst a + 1 .. fst b - 1|] |> Array.map (fun i -> (i, snd a))
            //let start = a + dirTemp

            let rList = rectList // Get a immutable snapshot of the current state
            // Find the index of the rectangle that contains both a and b
            let p = (-1, pointIsInRectangle rList a) ||> Array.fold (fun acc elem ->
                if (pointIsInRectangle rList b |> Array.fold (fun acc2 elem2 -> acc2 || (elem2 = elem)) false) then
                    elem
                else
                    acc)
            let splitRect = SplitPathBetweenPoints a b rectList.[p]
            rectList <- Array.append rectList [| (Array.append (fst splitRect) chordToAdd) |]
            rectList <- Array.append rectList [| (Array.append (snd splitRect) (Array.rev chordToAdd)) |]
            //printfn "%A" rectList
            rectList <- rectList.[p + 1 .. rectList.Length - 1]

        printf "chordless rectangulate"
        let jimbob = rectList |> Array.map (fun chordlessRect ->
            ChordlessRectangulate chordlessRect tiledShape) |> Array.collect (fun x -> x)
        printf "more collecting"
        jimbob
        